/* ClassicalGraphDefinition.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined(_ClassicalGraphDefinition_)
#define _ClassicalGraphDefinition_
#include <src/miscellaneous/OS_basic.h>
#include <src/data/Operators.h>
#include <src/algorithm/LinearAlgebra/VanDerMonde.h>
#include <src/data/TUple.h>
/*!
	OTSCentered structure tag
*/
struct OTSCentered{};
/*!
	Forward structure tag
*/
struct Forward{};
/*!
	Backward structure tag
*/
struct Backward{};
/*! Common Base class for graph definitions.
	Order is the Order of relations : i.e. the graph will be
	- 
	- Type = Forward	case : \f$ \nu_i = \{i-Order+1,...,i+Order\}\f$
	- Type = OTSCentered	case : \f$ \nu_i = \{i-Order,...,i+Order\}\f$
	- Type = Backward	case : \f$ \nu_i = \{i-Order,...,i+Order-1\}\f$
*/
template <class I
		, std::size_t Order 
		, class Centered_type 
		, class Boundary 
>
struct operator_interface_base_common
{
	enum:std::size_t									{Order = Order};
	typedef Centered_type								Centered_type;
	typedef Boundary									Boundary_type;	
	typedef typename I::argument_base_type argument_base_type;
	typedef typename argument_base_type::size_type size_type;
	typedef typename argument_base_type::base_value_type base_value_type;
	typedef TUple<base_value_type,Order+1> signature_type;
	operator_interface_base_common(I & i, const signature_type & sig) : call_back_(i), size_(0), sig_(sig) {
		typedef typename I::state_variable_type state_variable_type;
		typedef typename state_variable_type::value_type value_type;
//		call_back_.set_state_variable(state_variable_type(new value_type));
// computing coeff
	};
	virtual void rebuild()
	{};
	virtual void update()
	{};
	virtual void set_size(size_type size) {
		if (size_ != size) {
			size_ = size;
			for (OS_int i = 0;i< sig_.size(); ++i) 
				beta_[i] = (base_value_type) (i- (base_value_type) Order/2.) / (base_value_type) size_;   // geometry : the unit cube
			beta_ = VanDerMonde<signature_type>(beta_,sig_);

			if (size_ > 0) {
				rebuild();
				call_back_.build_graph();
			};
		};
	};
	size_type size() const {return size_;};

protected :
	size_type size_;
	I & call_back_;
	signature_type sig_;
	signature_type beta_;
};

template <class I, std::size_t O, class C, class B, std::size_t Dim = I::state_variable_type::value_type::Dim>
struct operator_interface_base : public operator_interface_base_common<I,O,C,B>
{
	OS_STATIC_CHECK(false); // to implement for your own set of data
};

template <class I, std::size_t O, class C, class B>
struct operator_interface_base <I,O,C,B,1> : public operator_interface_base_common<I,O,C,B>
{
	typedef operator_interface_base_common<I,O,C,B> Base;
	operator_interface_base(I & i, const signature_type & sig) : Base(i,sig){};
	virtual void update() {
		set_size(call_back_.get_state_variable()->size());
	};
};

template <class I, std::size_t O, class C, class B>
struct operator_interface_base <I,O,C,B,2> : public operator_interface_base_common<I,O,C,B>
{
	typedef operator_interface_base_common<I,O,C,B> Base;
	operator_interface_base(I & i, const signature_type & sig) : Base(i,sig){};
	virtual void rebuild() {
		OS_DYNAMIC_CHECK(false,"to implement in derived class");
	};
	virtual void update() {
		I::state_variable_type s_v = call_back_.get_state_variable();
		if (s_v->size(1) != size_) {
			size_ = s_v->size(1);
			for (OS_int i = 0;i< sig_.size(); ++i) 
				beta_[i] = (base_value_type) (i- (base_value_type) Order/2.) / (base_value_type) size_;   // geometry : the unit cube
			beta_ = VanDerMonde<signature_type>(beta_,sig_);

			if (size_ > 0)
				rebuild();
				call_back_.build_graph();
		};
	};

};
template <class I, std::size_t O, class C, class B>
struct operator_interface_base <I,O,C,B,3> : public operator_interface_base<I,O,C,B,2>
{
	typedef operator_interface_base<I,O,C,B,2> Base;
	operator_interface_base(I & i, const signature_type & sig) : Base(i,sig){};
	void rebuild() 
	{};
};
/*! Common non compilable Base class for graph definitions.
*/
template <class A
		, std::size_t Order 
		, class OTSCentered
		, class Boundary 
>
struct operator_interface : public operator_interface_base<A, Order, OTSCentered ,Boundary>
{
	OS_STATIC_CHECK(false);
//! Compile time error means that you have to implement your set of data
	operator_interface(){};
};

template <class I, std::size_t Order>
struct operator_interface<I,Order,OTSCentered,Periodic>
	: public operator_interface_base<I,Order,OTSCentered,Periodic>
{
	typedef operator_interface_base<I,Order,OTSCentered,Periodic> Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
	template <class T>
	void construct(T&) {
		OS_STATIC_CHECK(false);//("you have to implement the rebuild function for your own weight set type");
	};
	template <>
	void construct(weight_set_default_def<typename I::base_value_type> & weight_set) {
		typedef typename I::indice_type indice_type;
		typedef typename I::argument_base_type data;
		weight_set.resize(size_,size_,false);

		for (OS_size i = 0 ; i < size_ ; ++i)
			for (OS_size j = 0; j<= Order ; ++j) {
				indice_type k = jmmath::mod<data::indice_type,data::size_type>(i+j-Order/2,size_);
				if ( k >=0 && k < size_ && beta_[j]!=0 )
					weight_set(i,k) = beta_[j];
			};
	};
	void rebuild() 
	{
		construct(call_back_.weight_set_);
	};
};

template <class I
		, std::size_t Order 
>
struct operator_interface<I,Order,Backward,Periodic>
	: public operator_interface_base<I,Order,Backward,Periodic>
{
	typedef operator_interface_base<I,Order,Backward,Periodic> Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
	template <class T>
	void construct(T&) {
		OS_STATIC_CHECK(false);//("you have to implement the rebuild function for your own weight set type");
	};
	template <>
	void construct(weight_set_default_def<typename I::base_value_type> & weight_set) {
		typedef typename I::indice_type indice_type;
		typedef typename I::argument_base_type data;
		weight_set.resize(size_,size_,false);
		for (OS_size i = 0 ; i < size_ ; ++i)
			for (OS_size j = 0; j<= Order ; ++j) {
				indice_type k = jmmath::mod<data::indice_type,data::size_type>(i+j-Order/2,size_);
				if ( k >=0 && k < size_ && beta_[j]!=0 )
					weight_set(i,k) = beta_[j];
			};
	};
	void rebuild() 
	{
		construct(call_back_.weight_set_);
	};
};
template <class I
		, std::size_t Order 
>
struct operator_interface<I,Order,Forward,Periodic>
	: public operator_interface_base<I,Order,Forward,Periodic>
{
	typedef operator_interface_base<I,Order,Forward,Periodic> Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
	template <class T>
	void construct(T&) {
		OS_STATIC_CHECK(false);//("you have to implement the rebuild function for your own weight set type");
	};
	template <>
	void construct(weight_set_default_def<typename I::base_value_type> & weight_set) {
		typedef typename I::indice_type indice_type;
		typedef typename I::argument_base_type data;
		weight_set.resize(size_,size_,false);
		for (OS_size i = 0 ; i < size_ ; ++i)
			for (OS_size j = 0; j<= Order ; ++j) {
				indice_type k = jmmath::mod<data::indice_type,data::size_type>(i+j-Order/2,size_);
				if ( k >=0 && k < size_ && beta_[j]!=0 )
					weight_set(i,k) = beta_[j];
			};
	};
	void rebuild() 
	{
		construct(call_back_.weight_set_);
	};
};

/*!
	\brief Specialization for one dimensional Dirichlet OTSCentered operators
*/
template <class I
		, std::size_t Order 
>
struct operator_interface<I,Order,OTSCentered,Dirichlet>
	: public operator_interface_base<I,Order,OTSCentered,Dirichlet>
{
	typedef operator_interface_base<I,Order,OTSCentered,Dirichlet> Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {
	};
	template <class T>
	void construct(T&) {
		OS_STATIC_CHECK(false);//("you have to implement the rebuild function for your own weight set type");
	};
	template <>
	void construct(weight_set_default_def<typename I::base_value_type> & weight_set) {
		typedef typename I::indice_type indice_type;
		typedef typename I::argument_base_type data;
		weight_set.resize(size_,size_,false);
		for (OS_size i = 0 ; i < size_ ; ++i)
			for (OS_size j = 0; j<= Order ; ++j) {
				indice_type k = (indice_type) (i-Order/2+j);
				if ( k >=0 && k < size_ && std::abs(beta_[j]) > THRESHOLD )  //needed for singular operator
					weight_set(i,k) = beta_[j];
			};
	};
	void rebuild() 
	{
		construct(call_back_.weight_set_);
	};
	
};



template <class I
		, std::size_t Order 
>
struct operator_interface<I,Order,Backward,Dirichlet>
	: public operator_interface_base<I,Order,Backward,Dirichlet>
{
	typedef operator_interface_base<I,Order,Backward,Dirichlet> Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
	template <class T>
	void construct(T&) {
		OS_STATIC_CHECK(false);//("you have to implement the rebuild function for your own weight set type");
	};
	template <>
	void construct(weight_set_default_def<typename I::base_value_type> & weight_set) {
		typedef typename I::indice_type indice_type;
		typedef typename I::argument_base_type data;
		weight_set.resize(size_,size_,false);
		for (OS_size i = 0 ; i < size_ ; ++i)
			for (OS_size j = 0; j<= Order ; ++j) {
				indice_type k = (indice_type) (i-Order/2+j);
				if ( k >=0 && k < size_ && beta_[j]!=0 )
					weight_set(i,k) = beta_[j];
			};
	};
	void rebuild() 
	{
		construct(call_back_.weight_set_);
	};
};
template <class I
		, std::size_t Order 
>
struct operator_interface<I,Order,Forward,Dirichlet>
	: public operator_interface_base<I,Order,Forward,Dirichlet>
{
	typedef operator_interface_base<I,Order,Forward,Dirichlet> Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
	template <class T>
	void construct(T&) {
		OS_STATIC_CHECK(false);//("you have to implement the rebuild function for your own weight set type");
	};
	template <>
	void construct(weight_set_default_def<typename I::base_value_type> & weight_set) {
		typedef typename I::indice_type indice_type;
		typedef typename I::argument_base_type data;
		weight_set.resize(size_,size_,false);
		for (OS_size i = 0 ; i < size_ ; ++i)
			for (OS_size j = 0; j<= Order ; ++j) {
				indice_type k = (indice_type) (i-Order/2+j);
				if ( k >=0 && k < size_ && beta_[j]!=0 )
					weight_set(i,k) = beta_[j];
			};
	};
	void rebuild() 
	{
		construct(call_back_.weight_set_);
	};
};


template <class I, std::size_t O>
struct operator_interface<I,O,OTSCentered,Unbounded>
	: public operator_interface<I,O,OTSCentered,Dirichlet>
{
	typedef operator_interface<I,O,OTSCentered,Dirichlet>		Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
};
template <class I, std::size_t O>
struct operator_interface<I,O,Backward,Unbounded>
	: public operator_interface<I,O,Backward,Dirichlet>
{
	typedef operator_interface<I,O,Backward,Dirichlet>		Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
};
template <class I, std::size_t O>
struct operator_interface<I,O,Forward,Unbounded>
	: public operator_interface<I,O,Forward,Dirichlet>
{
	typedef operator_interface<I,O,Forward,Dirichlet>		Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
};

template <class I, std::size_t O>
struct operator_interface<I,O,OTSCentered,Neumann>
	: public operator_interface<I,O,OTSCentered,Dirichlet>
{
	typedef operator_interface<I,O,OTSCentered,Dirichlet>		Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
};
template <class I, std::size_t O>
struct operator_interface<I,O,Backward,Neumann>
	: public operator_interface<I,O,Backward,Dirichlet>
{
	typedef operator_interface<I,O,Backward,Dirichlet>		Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
};
template <class I, std::size_t O>
struct operator_interface<I,O,Forward,Neumann>
	: public operator_interface<I,O,Forward,Dirichlet>
{
	typedef operator_interface<I,O,Forward,Dirichlet>		Base;
	operator_interface(I & i, const signature_type & sig) : Base(i,sig) {};
};
#endif
